#---
#SIMULATION STUDY USING PARALLEL PROGRAMMING
#---

library(MASS)  

#run slaves
library(Rmpi)
mpi.spawn.Rslaves()

#compute MISE for several smoothing methods
par.sim.mise <- function(j, nid, m, tuning, T, s, r, bwSeq)
{
	source("sim_poisson_process.r")
	source("int_ems.r")
	source("int_local_constant.r")
	source("int_local_linear.r")

	#weights
	wgt <- c(0.066671344308688,
		0.149451349150581,
		0.219086362515982,
		0.269266719309996,
		0.295524224714753,
		0.295524224714753,
		0.269266719309996,
		0.219086362515982,
		0.149451349150581,
		0.066671344308688)

	#true intensity function
	int.true <- function(mesh, m, T, s, r)
	{
		m * dgamma(mesh, shape = s, rate = r)/pgamma(T, shape = s, rate = r)
	}

	#simulate data
	repeat
	{
		sim.data.temp <- sim.nhpp.data(nid, m, tuning, T, s, r)
		em <- ems.nhpp(panel.data = sim.data.temp$int, bw = 0, quad.order = 10)
		Jlength <- diff(sort(union(sim.data.temp$int$L, sim.data.temp$int$R)))
		lc.temp <- try(local.constant.nhpp(panel.data = sim.data.temp$int, bw = bwSeq[1], quad.order = 10), TRUE)
		ll.temp <- try(local.linear.nhpp(panel.data = sim.data.temp$int, bw = bwSeq[1], quad.order = 10), TRUE)

		if((length(wgt) * length(Jlength) == length(em$lattice)) & !is.character(lc.temp) & !is.character(ll.temp))
		{
			sim.int.data.lst <- sim.data.temp$int
			sim.exact.data.lst <- sim.data.temp$exact
			break
		}
	}		

	#compute MISE
	mise.exact <- rep(NA, length(bwSeq))
	mise.em_sl <- rep(NA, length(bwSeq))
	mise.em_sm <- rep(NA, length(bwSeq))
	mise.em_sr <- rep(NA, length(bwSeq))
	mise.lc <- rep(NA, length(bwSeq))
	mise.ll <- rep(NA, length(bwSeq))

	#exact, left, mid and right points
	em <- ems.nhpp(panel.data = sim.int.data.lst, bw = 0, quad.order = 10)
	midpts <- sort(union(sim.int.data.lst$L, sim.int.data.lst$R))
	rendpts <- midpts[-1] 
	lendpts <- midpts[-length(midpts)]
	Jlength <- diff(midpts)
	midpts <- (lendpts + rendpts)/2

	for(k in 1:length(bwSeq))
	{
		#exact
		exact <- density(sim.exact.data.lst$TIME, bw = bwSeq[k], from = 0, to = T)
		mise.exact[k] <- sum((length(sim.exact.data.lst$TIME)/nid * exact$y - int.true(exact$x, m, T, s, r))^2 * diff(exact$x[1:2]))
    
		#EM + S (left)
		em_s <- as.vector(em$Lambda.est %*% outer(lendpts, em$lattice, function(x, y) dnorm(x, mean = y, sd = bwSeq[k])))
		mise.em_sl[k] <- as.vector(as.vector(outer(wgt, Jlength/2)) %*% (em_s - int.true(em$lattice, m, T, s, r))^2)

		#EM + S (mid)
		em_s <- as.vector(em$Lambda.est %*% outer(midpts, em$lattice, function(x, y) dnorm(x, mean = y, sd = bwSeq[k])))
		mise.em_sm[k] <- as.vector(as.vector(outer(wgt, Jlength/2)) %*% (em_s - int.true(em$lattice, m, T, s, r))^2)
    
		#EM + S (right)
		em_s <- as.vector(em$Lambda.est %*% outer(rendpts, em$lattice, function(x, y) dnorm(x, mean = y, sd = bwSeq[k])))
		mise.em_sr[k] <- as.vector(as.vector(outer(wgt, Jlength/2)) %*% (em_s - int.true(em$lattice, m, T, s, r))^2)
     
		#local constant local EM
		lc <- local.constant.nhpp(panel.data = sim.int.data.lst, bw = bwSeq[k], quad.order = 10)
		mise.lc[k] <- as.vector(as.vector(outer(wgt, Jlength/2)) %*% (lc$lambda.est - int.true(lc$lattice, m, T, s, r))^2)
    
		#local linear local EM
		ll <- local.linear.nhpp(panel.data = sim.int.data.lst, bw = bwSeq[k], quad.order = 10)
		mise.ll[k] <- as.vector(as.vector(outer(wgt, Jlength/2)) %*% (ll$lambda.est - int.true(ll$lattice, m, T, s, r))^2)
	}

	return(list(sim.int.data.lst, sim.exact.data.lst, 
		mise.exact, mise.em_sl, mise.em_sm, mise.em_sr, mise.lc, mise.ll))
}

#true shape and rate parameters
s <- 4.75
r <- 0.75

#set.seed(100)
#simulate exact and censored times for panel count data
nid <- 30 #number of subjects per dataset
m <- 5 #number of event times per subject
T <- 20 #number of visiting times
nDataset <- 500 #number of simulated datasets
tuning <- 0.1 #parameter to adjust probabilities of missing a visit (between 0 and 1)

#bandwidth
bwSeq <- seq(0.05, 2.45, by = 0.05)

#compute data and MISE
mise.temp <- mpi.parSapply(1:nDataset, par.sim.mise, nid, m, tuning, T, s, r, bwSeq, simplify = FALSE)

save.image("sim_pc30_100.RData")

sim.int.data.lst <- NULL
sim.exact.data.lst <- NULL
mise.exact <- NULL
mise.em_sl <- NULL
mise.em_sm <- NULL
mise.em_sr <- NULL
mise.lc <- NULL
mise.ll <- NULL
err <- c()
for(j in 1:length(mise.temp))
{
	if(is.character(mise.temp[[j]]))
	{
		err <- c(err,j)
	}
	
	else
	{
		sim.int.data.lst <- c(sim.int.data.lst, list(mise.temp[[j]][[1]]))
		sim.exact.data.lst <- c(sim.exact.data.lst, list(mise.temp[[j]][[2]]))
		mise.exact <- rbind(mise.exact, mise.temp[[j]][[3]])
		mise.em_sl <- rbind(mise.em_sl, mise.temp[[j]][[4]])
		mise.em_sm <- rbind(mise.em_sm, mise.temp[[j]][[5]])
		mise.em_sr <- rbind(mise.em_sr, mise.temp[[j]][[6]])
		mise.lc <- rbind(mise.lc, mise.temp[[j]][[7]])
		mise.ll <- rbind(mise.ll, mise.temp[[j]][[8]])
	}
}

#compute average MISE
avg.mise.exact <- apply(mise.exact, 2, mean)
avg.mise.em_sl <- apply(mise.em_sl, 2, mean)
avg.mise.em_sm <- apply(mise.em_sm, 2, mean)
avg.mise.em_sr <- apply(mise.em_sr, 2, mean)
avg.mise.lc <- apply(mise.lc, 2, mean)
avg.mise.ll <- apply(mise.ll, 2, mean)

#computer average probability of missing a visit between 0 and T
avg.miss.prob <- mean(((1:T)/T)^tuning - 0.05)

save.image("sim_pc30_100.RData")
mpi.close.Rslaves()
mpi.quit()
